Random Sampling vs. Exact Enumeration of Attractors in Random Boolean Networks 
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We clarify the effect different sampling methods and weighting schemes have on the statistics of 
attractors in ensembles of random Boolean networks (RBNs). We directly measure cycle lengths 
of attractors and sizes of basins of attraction in RBNs using exact enumeration of the state space. 
In general, the distribution of attractor lengths differs markedly from that obtained by randomly 
choosing an initial state and following the dynamics to reach an attractor. Our results indicate that 
the former distribution decays as a power-law with exponent 1 for all connectivities A" > 1 in the 
infinite system size limit. In contrast, the latter distribution decays as a power law only for K = 2. 
This is because the mean basin size grows linearly with the attractor cycle length for K > 2, and is 
statistically independent of the cycle length for K = 2. We also find that the histograms of basin 
sizes are strongly peaked at integer multiples of powers of two for K <?>. 

PACS numbers: 02.70.Uu, 05.10.Ln, 87.10.+e, 89.75.Fb, 89. 75. He 
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I. INTRODUCTION 

Random Boolean Networks (RBNs) [l| have been 
widely used as elementary models for genetic regula- 
tion. In such a model, a binary state based on Boolean 
logic encapsulates local gene expression. An RBN 

consists of N Boolean (0,1) elements where the value of 
each element evolves in discrete time according to a ran- 
dom Boolean function of K randomly chosen distinct in- 
puts. Within an annealed approximation [5], the RBN 
can be in one of three phases: a frozen phase {K = 1), 
in which a perturbation to the state of a single node can 
propagate only to a finite number of nodes; a chaotic 
phase {K > 2), in which the perturbation spreads to a 
finite fraction of the nodes; and a critical phase {K — 2), 
which lies in between the frozen and chaotic phases. 

In a finite RBN the number of possible states is also 
finite. Therefore, each dynamical state in a determin- 
istically updated RBN is either transient or belongs to 
an attractor cycle. A transient state may be reached no 
more than once during the dynamics. Meanwhile, a state 
which belongs to an attractor cycle may be reached in- 
finitely often if the chosen initial state falls into the basin 
of attraction of the attractor cycle. Many analytical re- 
sults concerning the duration of transients and attractor 
cycle lengths have been established for the random map 
(RM), which is the limit of RBNs when K N. These 
analytical arguments are based on the fact that the an- 
nealed approximation is exact for the random map 0, 0] ■ 
However, most of the known results for K = 2 critical 
RBNs come from numerical simulations. 

Since the size of the state space grows as 2^, an ex- 
haustive search of attractors is infeasible for large N. 
Instead, a random sampling procedure has (almost ex- 
clusively) been used jj, ^, 3] ■ the system is prepared in a 
randomly chosen initial state, and the RBN rules are used 
to evolve the system until an attractor is reached. Typi- 
cally a fixed number of initial conditions are used for each 
realization of an RBN, and the ensemble properties are 
obtained by repeating this procedure for many realiza- 



tions. Using this sampling procedure, it has been shown 
that the lengths of the transients to reach an attractor, 
as well as the lengths of the attractor cycles reached are 
both power-law distributed for the ensemble of K = 2 
critical RBNs Q. However, the size of the state space 
limits the accuracy of the random sampling procedure. 
For instance, small basins of attraction may be under- 
sampled as suggested by recent studies of the (mean) 
number of attractors [13, [HI [l3| ■ This has raised strong 
concerns regarding the validity of estimates of the distri- 
butions of cycle lengths and sizes of basins of attraction 
that are based on the random sampling procedure [J — 
like the ones obtained in Ref. 

In particular, the simple sampling procedure described 
above can only measure the distribution of attractor 
lengths reached from a randomly chosen initial state, but 
not the unbiased distribution of attractor cycle lengths 
itself. These two distributions only coincide if there are 
no correlations between the length of an attractor cycle 
and its basin size. While it can be shown analytically 
that the unbiased distribution of attractor cycle lengths 
decays as a power-law with exponent 1 for the random 
map Q, less is known about K < N. Refs. [isl [l^ esti- 
mate a basin entropy using exact enumeration for various 
K , but do not consider attractor cycle lengths. 

In this article, we clarify the effects of different sam- 
pling procedures and weighting schemes by making exact 
enumerations of the state space of RBNs to estimate the 
distributions of both the sizes of the basins of attrac- 
tion and the attractor cycle lengths. By weighting each 
attractor by its basin size, we can reproduce the distri- 
butions obtained by randomly sampling initial states as 
discussed above. We find numerically that the unbiased 
distribution of attractor cycle lengths differs markedly 
from that obtained by randomly sampling initial states 
for all K > 2. This is corroborated by analytical argu- 
ments based on an annealed approximation [7] . Remark- 
ably, for K — 2 both distributions are well approximated 
by the same power-law. This difference between criti- 
cal and chaotic RBNs is related to how the basin size 
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depends on the length of its attractor cycle. We show 
analytically that the mean basin size increases linearly 
with the length of the attractor for the random map. 
Numerically, we find that this also holds for K > 2. For 
K — 2, the mean basin size is independent of the attrac- 
tor cycle length. Enumeration also allows us to study 
the distribution of basin sizes; we find that for K = 2, 
those distributions are discrete and strongly peaked at 
integer multiples of powers of two, reflecting in part the 
analytical structure previously found ior K = 1 \W^ . 

In Section II, we present a physical motivation for the 
different weighting schemes used to compute various dis- 
tributions for cycle lengths or basin sizes as well as a 
more formal discussion clarifying the mathematical rela- 
tions between these distributions. Section III contains 
the results from numerical simulations, while Section IV 
concludes with a summary of the main findings. 



the notation we use, Q denotes distributions obtained by 
weighting attractors according to the size of their basin of 
attraction, while P denotes distributions obtained with- 
out regard to the basin size. In addition, the subscript 
u in Qu{l) denotes a distribution obtained by weighting 
each RBN equally (e.g. making a histogram for each 
RBN and then uniformly averaging the results over dif- 
ferent realizations), while the subscript a in Pa{l) denotes 
a distribution in this case obtained by weighting each at- 
tractor equally - regardless of which RBN it came from 
(see below for further discussion). It is obvious that there 
are many other distributions that can be estimated, de- 
pending on the weighting scheme. In general, there is 
no reason these distributions should be similar. How- 
ever, there are mathematical relations connecting them 
as described next. 



II. DEFINITION OF DISTRIBUTIONS 

Depending on the weighting scheme, different distribu- 
tions for the same quantity can be obtained. One can, for 
instance, make an estimate of the distribution of attrac- 
tor lengths that would be obtained by randomly sampling 
/ initial states (for each of R realizations of the RBN) and 
following the dynamics to reach an attractor as follows: 
For an RBN, make an exhaustive list of each attractor 
with its cycle length and size of its basin. Then pick / 
attractors randomly - each with a weight proportional 
to its basin size - and compute a histogram of attractor 
cycle lengths for the RBN. Repeat this for R different re- 
alizations of the RBN and average the results to obtain 
an ensemble averaged probability distribution. 

The distribution of cycle lengths obtained this way cor- 
responds to an estimate of Qu{l) in the notation below. 
Note that using a weight proportional to the basin size 
accounts for the fact that initial states have a proportion- 
ally higher probability to fall in larger basins within the 
state space than in smaller ones. For simplicity, we will 
refer to Qu{l) as the distribution of cycle lengths obtained 
by random sampling - even though it oversimplifies the 
situation. In particular, as this discussion implies, we are 
not using a random sampling method, but rather repro- 
ducing what would be obtained using that method [23| . 

On the other hand, using the exhaustive list of each at- 
tractor and its basin size for all R realizations of an RBN, 
one could compute a histogram by directly accumulating 
the results for each attractor on the list - independent of 
its basin size and also independent of which RBN realiza- 
tion it appears in within the ensemble of R realizations. 
Such a distribution of attractor lengths corresponds to 
Pa{l) in the notation below. In order to assist the gen- 
eral reader, we will refer to Pa{l) as the distribution of cy- 
cle lengths obtained by exact enumeration ~ again, even 
though it oversimplifies the true situation. 

These two distributions differ in two ways as signified 
by the P vs. Q label as well as the subscript a vs. u. In 



A. Formal development 

The state space of a single realization of an RBN may 
contain many different attractors. Each attractor is char- 
acterized by {l,h), which are its cycle length / and the size 
h of its basin of attraction. For RBN « in a given ensem- 
ble, we count each attractor a in its state space, and 
record the respective / and h. This allows us to obtain 
the probability that a randomly chosen attractor of RBN 
i has cycle length / and basin size 6, 

. A, 

a— 1 

A, ' 

Here Ai{l, b) is the number of attractors with cycle length 
I and basin size b in RBN i, and Ai is the total number of 
attractors in it. The probability that a randomly chosen 
state lies in the basin of an attractor with {I, b) for RBN 
i is given by 

a=l ^ 

= ^A,[l,h) = ^MP^^{l,b). (2) 

It is easy to check that both P^''''{l,b) and Q'^^^l^b) are 
normalized to unity. 

To obtain the corresponding probabilities for attractor 
lengths and basin sizes for an ensemble of RBNs, we as- 
sign a normalized weight Wi to each RBN i. If i? is the 
number of RBNs in the ensemble, the probability that a 
randomly chosen attractor from the ensemble has cycle 
length / and basin size 6, is 

R 

P(Z,6;K}) = ^w,P«(Z,&). (3) 

i=l 
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Similarly, the probability over an ensemble of RBNs that 
a randomly chosen state lies in the basin of an attractor 
with b) is given by 



Qil,b;{w,}) = ^u;,QW(/,6) 



2N 



P{l,b;{w,A,}). 



(4) 



For uniform weights Wi = l/R, we obtain the distribu- 
tion, 



1 Y^ Mi^b) 



(5) 



If, in contrast, we use the weights Wi = Ai / {R{A)) , we 
obtain a distribution for the ensemble where all attractor 
are counted equally: 



(6) 



where (A) is the mean number of attractors in a single 
realization of an RBN. 

Obtaining ensemble averaged "Q" distributions pro- 
ceeds in the same manner. For uniform weights Wi — 
1/R, Eq. (g]) becomes 



= —{A)Pa{l,b) 



(7) 



where we used Eq. ([6]) to obtain the second equality. 

For comparisons, we also consider Eq. ^ when the 
RBNs are weighted by the inverse of the number of at- 
tractors, 

Qi/a{l,b) EE Q{l,b;{A-'/{R{A-'))}) 



1 ^ 
1 X - 

4-l\ W ^ 



Ail,b) 



i?(A-i)2^^ A, 
1 b 



(8) 



Eqs. dSEHfe El) make up the four weighting schemes we 
focus on in the remainder of this paper. 

Using these joint distributions, we can construct the 
distributions of attractor lengths, by summing over the 
basin sizes b. For example, 

Q„(Z) = ^Q„(Z,6) = ^^5P,G,6) 



2N 



{b{l))aPail). 



(9) 



Here we have used Eq. ^ and defined the mean basin 
size of attractors of length I, 



{b{l))a 



~ Pail) ' 



(10) 



with 



Pail) EE Y.Pail,b) 



(11) 



The distributions Puil) and Qi/a(0 ^''^ defined anal- 
ogously. In order to assist the general reader, we refer 
to Quil) as the distribution of cycle lengths obtained by 
random sampling and Pail) as the distribution obtained 
by exact enumeration. 

In this paper, all these quantities are estimated by ex- 
act enumeration of the state space. For each realization 
of the RBN, we find the image of each of the 2^ states 
under the dynamical map. We connect each state and 
its image by a directed link to form the state space net- 
work (SSN). We follow these directed links to reach the 
attractors, which are cycles of directed links. For each 
attractor we find its cycle length and basin size p^ . 



III. RESULTS 

Fig. [1] shows the distribution of attractor lengths ob- 
tained by random sampling Quil), the distribution ob- 
tained by exact enumeration Pail), as well as Puil) and 
Qi/ail) ior K = 1, K = 2 and the random map. Re- 
sults for other values of K were also obtained and some 
of these results are shown later. We find that for all 
K ^ 1, Puil) is similar to Pail), while Quil) is similar to 
Qi/ail)- Defining AiH) = J2b b), these results imply 
that the ratio Aiil)/Ai is statistically independent of the 
total number of attractors in the RBN, Ai, for K ^ 1. 
This follows from comparing Eqs. (O and © or their 
analogs for Qu{l) and Qi/ail), respectively. Such inde- 
pendence is expected for the random map and for K > 2 
in the thermodynamic limit, based on the annealed ap- 
proximation We currently have no explanation for 
why this would also hold ioi K — 2. For K — 1, the dif- 
ference between Puil) and Pail) or Quil) and Qi/a(0 is 
a direct indication that there are significant correlations 
between Aiil)/Ai and Ai. This is expected since long 
attractors arise from long loops in the network of RBN 
elements, which also give rise to a larger total number of 
attractors [H, [l^ . Since we are mostly concerned with 
ii' > 1 in what follows, we will focus on the distribution 
of cycle lengths obtained by exact enumeration. Pail), 
and the distribution of cycle lengths obtained by ran- 
domly sampling, Quil), which differ from each other for 
all A' 7^ 2 as shown in Fig. [1] 

While this difference is well-known for the random 
map, the relationship between these two distributions has 
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FIG. 1: (Color online.) Distribution of attractor cycle lengths using four weighting schemes for K — 1, K = 2 and the 
random map. Results are statistically indistinguishable for K = 2, while differences appear for all other K. Qu{l) represents 
the (simplest) distribution of cycle lengths obtained by random sampling, while Pa(l) represents the (simplest) distribution 
obtained by exact enumeration - as explained in Section II. This and all subsequent figures are based on 5 x 10^ RBN realizations 
with N = 22, unless stated otherwise. 
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not been clarified for other values of K. As shown for the 
random map in Refs. [gI. [T7|. 



10 - 



Quil) 
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(2^ - 1)! 



{2N -t)l{2^y 

2N/2 
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while in Refs. [1, [iBl it was found that 



Pail) 



(12) 



(13) 



The latter distribution decays as a power-law with expo- 
nent 1, up to an dependent cut-off, while the former 
- being the complement of the error function - is a flat 
distribution up to an N dependent cut-off. Indeed, us- 
ing the finite size scaling method for different N allows 
one to collapse the different curves for Pa{l) and Qu{l), 
respectively. As discussed in Section II, Q„(/) is related 
to Pa{l) by a factor which is the mean basin size of at- 
tractors of length I, (&(/)) a, see Eq. (0). Hence Qu(0 ^^'^ 
Pa{l) can only have different functional forms if {b{l))a 
varies with I. 

From Eqs. (fT^ and p^ . it is straightforward to show 
that {b{l))a / for the random map. Fig. [2] confirms this. 
Moreover, Fig. [5] shows that the linear relationship also 
holds for large K suggesting that the qualitative differ- 
ences between Pa{l) and Qu{l) persist. Indeed, using the 
annealed approximation, Eqs. ((T2)) and (fT3)) can be ex- 
tended to chaotic RBNs {K > 2), as was done previously 
for Eq. in [7] . This leads to a slightly modified finite 
size scaling form. For AT > 6, the expressions given in 
Ref. allow us to collapse the different curves as shown 
in Fig. [31 For smaller K the accessible values of A^ are 
too small to apply the annealed approximation, which is 
exact in the limit N ^ oo. This is also confirmed by 
Fig. O which indicates that the dependence of {b{l))a on 
I weakens with decreasing K for fixed A^. This explains 
why for N = 22 and AT = 3, Pa{l) and Qu{l) do not show 
the same qualitative differences as seen for higher values 
of a: (see Fig. 111). 

Fig. [T] also indicates that Pa{l) and Qu{l) are statisti- 
cally identical for the critical case K = 2. As discussed in 
Section II (see Eq. (O), this can be directly related to the 
observation in Fig.[2]that {b{l))a is basically independent 
of Z for AT = 2. Numerical results presented in Ref. 
suggest that Qui^) for AT = 2 decays as a power-law with 
an A^-dependent exponent that approaches 1 in the ther- 
modynamic limit. Our results in Fig. 0] show that the 
exponent of Pa{l) also varies with system size. Provided 
that {b{l))a remains independent of I for N oo, this 
implies that the exponent for Pa{l) also approaches 1 in 
this limit. 

Consequently, the difference between the state spaces 
of chaotic and critical RBNs lies not in the distribution 
of cycle lengths obtained by exact enumeration, Pa{l) 
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FIG. 2: (Color online.) Mean (normalized) basin size, 
{6(Z))a/2^, vs. attractor length, /, for various values of K and 
the random map (RM). {b{l))a decreases with I for A = 1, 
and increases for K > 2 and the random map. For K — 2, 
the mean basin size is approximately independent of I. Note 
that the estimated mean basin size is only shown for those 
bins which had at least 50 data points. 



(Fig. [5l upper panel), which decays as a power law for 
all A" > 1, but in the correlations between the attractor 
cycle length and its basin size. While {b{l))a ~ constant 
for ensembles of AT = 2 critical RBNs, {b{l))a oc I for 
ensembles of chaotic RBNs. The presence or absence of 
such correlations is naturally captured by whether or not 
Qu{l) (Fig.El lower panel) differs from Pa{l)- 

Most approaches to estimate attractor length distri- 
butions in RBNs in the past have exclusively focused 
on Qu{l), which is the simplest distribution that can 
be obtained by randomly sampling initial states. In- 
deed, the specific random sampling methods applied in 
Refs. [ll, simply not able to measure {b{l))a or 

Pa{l)- This is the clear advantage of the exact enumera- 
tion method applied here. 

Exact enumeration also allows us to find the sizes of 
the basins of attraction. The distribution of basin sizes is 



known analytically for K = 1 .15.1 . AT = A^ [6| and all K 
in the chaotic phase • Ref. [14| used exact enumeration 
to numerically find that the basin entropy scales only for 
AT = 2 critical RBNs. However, they did not present the 
distribution of basin sizes for AT = 2. 

Fig. [5] shows the distribution of basin sizes Pa{b) = 
J2i Pa{l, b) for a: = 1 and AT = 2. For AT = 1 the basin 
size is given by 
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(14) 



where m is the number of relevant elements in the 
RBN [l5|. This causes the distribution of basin sizes 
for A' = 1 to be discrete and peaked at integer multiples 
of powers of two. In the lower panel of Fig. [S] we see that 
the distribution for A' = 2 shares some of this special 
structure - it is dominated by, but not entirely restricted 



6 



10' 



10' 



.10" 



CL _2 

10 



10" 







^N=10 






■e-N=l3 






-i-N=16 






-B-N=19 






■A'N=22 


10° 




^t*^^^ >N=25 




10-^ 






^«io-* 






- 10-° 






10° 10' io" io' lo" 

1 





10" 



10 10 

|/((2N)0.42, 



10" 



10' 




(a) 




10"^ 10"^ 10"^ 10° 10^ 

|/(2N)0.42 

(b) 

FIG. 3: (Color online.) Finite size scaling analysis of attractor 
length distributions for /sT = 6 RBNs, (a) Pa{l) and (b) Qu{l). 
The insets show the original distributions. These results are 
based on the following number of realizations: 5 x 10'' for 
iV = 10 - 22, 5 X 10'' for iV = 25 and 5 x 10^ for N = 28. 

to, integer multiples of powers of 2. This is most likely 
because the structure of if = 2 relevant components is 
similar to those of if = 1 as discussed in [2]. Approxi- 
mating Pa{b) for K = 2 hy & continuous distribution is 
shown in the inset of Fig. [51 

IV. CONCLUSION 

The full enumeration of the state space for each real- 
ization of an RBN presented here allows us to estimate 
the distribution of attractor cycle lengths or basin sizes 
mimicking different sampling procedures and weighting 
schemes. The unbiased distribution of attractor lengths, 
obtained by weighting all attractors in the ensemble 



FIG. 4: (Color online.) Pa{l) for critical K ^2 RBNs with 
various A'^. The slope of the curves varies systematically with 
A^. Based on the results of @] we expect this distribution to 
be a power law with exponent —1 in the thermodynamic limit. 
These results are based on 5 x 10^ realizations for A'^ — 10 — 22, 
5 X 10'' for iV = 25 and 5 X 10^ for = 28. 



equally, is very different from the cycle length distribu- 
tion of attractors reached from a randomly chosen initial 
state. Yet, for the critical case K = 2 both distribu- 
tions are statistical identical. This directly implies that 
the random sampling procedure applied in Ref. 0] for 
K — 2 and large N indeed allows one to obtain a reli- 
able estimate for the unbiased distribution of attractor 
cycle lengths. The comparison of the different weight- 
ing schemes also shows that the fraction of attractors of 
a given length in a given RBN is statistically indepen- 
dent of the total number of attractors in that RBN for 
all K >1. 

In addition, our findings show that the existence of a 
power-law decay in the unbiased distribution of attractor 
lengths is not an indicator of the criticality of an RBN 
(as defined in the introduction), since this distribution 
also decays as a power law for all RBNs in the chaotic 
phase {K > 2). Thus, the difference between the state 
spaces of chaotic and critical RBNs lies not in the un- 
biased distribution of attractor cycle lengths, but in the 
correlations between the attractor cycle length and the 
sizes the basins that these attractors drain. The typi- 
cally applied random sampling procedure naturally cap- 
tures this correlation between the attractor length and 
the basin size. This correlation is, however, only made 
explicit by measuring the joint probability of attractor 
length and basin size, which is accomplished by a full 
enumeration of the state space of an RBN. Such an enu- 
meration scheme allows one in particular to obtain the 
distribution of basin sizes. For K=2, we find that the dis- 
tribution is strongly peaked at integer multiples of powers 
of 2. In a more general setting, this enumeration scheme 
could also prove useful to investigate correlations between 
attractor cycle length and basin size in systems like cellu- 
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lar automata [T^ , discrete dynamical mappings [20j 
and multi-stable dynamical systems l2lll22l. 23 | . or in 
models of genetic regulatory networks [2J, [2^ • 



10" 



010" 



10" 



10" 



10' 



(a) 







X X — X )( )<W)( x~ 

-eK=i \ 

*K=2 \ 
-B-K=3 
-*-K=4 
-0-K=6 





10' 

I 



lO'' 



10' 



(b) 



ACKNOWLEDGMENTS 



FIG. 5; (Color online.) (a) Pa{l) and (b) Qu{l) for various 
K and the random map (RM). Pa{l) is a power law for all 
K > 1. In the thermodynamic limit all > 1 distributions 
are expected to have the same exponent, I. Qu{l) is a power 
law only for the critical K = 2 distribution. Note that the 
curve for the critical K — 2 ensemble does not exhibit an 
apparent cut-off in either case. 



We thank P. Grassberger for helpful discussions and 
useful comments. This work was partially supported by 
NSERC. 



[1] 
[2] 



[3] 
[4] 



[5] 
[6] 
[7] 



S. Kauffman, J. Theor. Biol. 22, 437 (1969). 
B. Drossel, Reviews of Nonlinear Dynamics and Com- 
plexity (Wiley, 2008), vol. 1, pp. 69-110, ISBN 978-3- 
527-40729-3. 

S. Bornholdt, Science 310, 449 (2005). 

M. Aldana, S. Coppersmith, and L. Kadanoff, 

In Perspectives and Problems in Nonlinear Science 

(Springer,New York, 2003), pp. 23-89. 

B. Derrida and Y. Pomeau, Europhys. Lett. 1, 45 (1986). 

B. Derrida and H. Flyvbjerg, J. Physique 48, 971 (1987). 

U. BastoUa and G. Parisi, Physica D 98, 1 (1996). 



[8] A. Bhattacharjya and S. Liang, Phys. Rev. Lett. 77, 1644 
(1996). 

[9] M. Paczuski, K. Bassler, and A. Corral, Phys. Rev. Lett. 
84, 3185 (2000). 
[10] B. Samuelsson and C. Troein, Phys. Rev. Lett. 90, 98701 
(2003). 

[11] S. Bilke and F. Sjunnesson, Phys. Rev. E 65, 016129 
(2001). 

[12] B. Drossel, T. Mihaljev, and F. Greil, Phys. Rev. Lett. 

94, 88701 (2005). 
[13] P. Krawitz and I. Shmulevich, Physical Review E 76, 



8 



10" 



CL 



0.4 
0.35 

0.3 
0.25 
\ 0.2 

L 

0.15 
0.1 
0.05 


0.06 
0.05 
0.04 
\ 0.03 

L 

0.02 
0.01 




10" 




K=1 



lO'' 



10' 



10^ lo'' 10^ 



11 12 13 14 15 16 17 18 19 20 21 22 
loggb 



K=2 




10^ 10^ 10^ lO'' lo'^ lo'' 



1 1 !■ I 



[14 
[15: 

[16: 

[17: 
[is: 

[19 

[20 

[21 

[22' 
[23: 

[24' 
[25 

[26 

[27 



36115 (2007). 

P. Krawitz and I. Shmulevich, Phys. Rev. Lett 98, 158701 

(2007) . 

H. Flyvbjerg and N. Kjaer, J. Phys. A: Math. Gen. 21, 
1695 (1988). 

A. Shreim, A. Berdahl, V. Sood, P. Grassberger, and 
M. Paczuski, New J. Phys. 10, 013028 (2008). 

B. Harris, Ann. Math. Statist 31, 1045 (1960). 
S. Wolfram, Rev. Mod. Phys. 55, 601 (1983). 

A. Shreim, P. Grassberger, W. Nadler, B. Samuelsson, 
J. Socolar, and M. Paczuski, Phys. Rev. Lett. 98 (2007). 
F. Kyriakopoulos and S. Thurner, Lecture Notes in Com- 
puter Science 4488, 625 (2007). 

U. Feudel, C. Grebogi, and B. R. Hunt, Physical Review 
E 54, 71 (1996). 

U. Feudel and C. Grebogi, Chaos 7, 597 (1997). 

U. Feudel and C. Grebogi, Phys. Rev. Lett. 91, 134102 

(2003). 

R. Albert and H. Othmer, J. Theor. Biol. 223, 1 (2003). 
A. Samal and S. Jain, BMC Systems Biology 2, 21 

(2008) . 

M. Davidich and S. Bornholdt, PLoS ONE 3, el672 
(2008). 

In principle, one could use a random sampling method 
to reproduce all the other distributions mentioned below 
by appropriate (re)weighting, e.g. with inverse basin size, 
etc., if the relevant quantities were measured. 
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FIG. 6: (Color online.) The unbinned (main plots) and 
binned (insets) distribution of basin sizes, Pa{b), for K = 1 &i 
2. For K = 1 the distribution is discrete with peaks only at 
integer multiples of powers of 2 as found in [l5| |. The distri- 
bution for K = 2 shares some of this special structure - it is 
dominated by, but not entirely restricted to, integer multiples 
of powers of 2. 



